






































AD-A144 216 


KFOSR-TR- 8 4-0 607 






AFOSR - m - 83 - 0189 


HARMONIC CONTROL TO REDUCE TORQUE PULSATIONS 
IN BRUSHLESS DC MOTOR DRIVES 


Jimmie J. Cathey 

University of Kentucky Research Foundation 
Lexington, KY 40506 


March 1984 


Final report for period May 1983 - January 1984 


8 


Approved for public release; Distribution unlimited 


LU 

fr 


f Prepared for 

AIR FORCE OFFICE OF SCIENTIFIC RESEARCH 
BOLLING AIR FORCE BASE, D.C. 20332 


s 


DTiC 

elected 

Aug 8 1984 


u 


84 08 06 129 



















TABLE OF CONTENTS 


SECTION PAGE 

I. INTRODUCTION 1 

1.0 Background 1 

1.1 Orientation 1 

1.2 Identification of Problem 2 

1.3 Consideration of Alternative Technologies 2 

1.3.1 Operation for Speeds Not Near Zero 2 

1.3.2 Operation for Speeds Near Zero 4 

1.3.3 Operation for Regenerative Power Flow 4 

2.0 Objective 5 

3.0 Scope 5 

II. PROCEDURE 7 

1.0 System Description 7 

2.0 Assumptions 7 

3.0 Mathematical Model 9 

4.0 Current Harmonics to be Eliminated 10 

5.0 Modulation of Phase Voltage 11 

5.1 Development of Modulation Function 11 

5.2 Minimization of Slack Variable 12 

5.3 Cycloconverter Connection Method 15 

5.4 Fourier Spectrum Analysis 15 

III. RESULTS 17 

1.0 Summary 17 

2.0 Performance Studies 18 

2.1 4500 RPM (10Z Speed Case) 19 

2.1.1 Unmodulated Phase Voltage 19 

2.1.2 Elimination of Sixth Harmonic of Torque 25 

2.1.3 Elimination of Sixth and Twelfth Harmonics of Torque 25 

2.2 2250 RPM (5% Speed Case) 31 

2.2.1 Unmodulated Phase Voltage 31 

2.2.2 Elimination of Sixth Harmonic of Torque 31 

2.2.3 Elimination of Sixth and Twelfth Harmonics of Torque 41 

2.2.4 Elimination of Sixth and Twelfth Harmonics of Torque 

(o - 15°) 41 

2.2.5 Elimination of Sixth and Twelfth Harmonics of Torque 

(a - 40°) 50 

2.3 450 RPM (1Z Speed Case) 59 

2.3.1 Unmodulated Phase Voltage 59 

2.3.2 Elimination of Sixth Harmonic of Torque 59 

2.3.3 Elimination of Sixth and Twelfth Harmonies of Torque 64 



.AIR F03Ti 
N0'TIC2'>? /- 
This tf-f 

epp-cTf.^ 
Distrit l 
■ATTEST J. Rl.:. 


Chief, Technical Information Division 








TABLE OF CONTENTS (CONTINUED) 


SECTION 

IV. DISCUSSION 

1.0 Increase In Ohmic Losses 
2.0 Cycloconverter Mode Change 

3.0 Extension to Higher Harmonics 


PAGE 

73 

73 

73 

74 


& 

I 


V. CONCLUSIONS AND RECOMMENDATIONS 


75 


VI. REFERENCES 


77 


APPENDICES 

A HARMONICS OF PWM TIME FUNCTION 79 

B MODULATION FUNCTION PROGRAMS 81 


1.0 

Initial Solution 

81 

1.1 

Theory of Harmonic Elimination 

81 

1.2 

Initial Solution Program 

84 

2.0 

Optimization of Modulation Function 

98 

2.1 

Discussion.of Procedure 

98 

2.2 

Modulation Function Optimization Program 

101 


C FOURIER SPECTRUM PROGRAM 117 


D PERFORMANCE PROGRAM 


123 




iv 








LIST OF ILLUSTRATIONS 


FIGURE PAGE 

1. Alternative Power Circuits 3 

2. Schematic of Power Circuit 8 

3. Approximate Form of Motor Phase Current 10 

4. Unmodulated Phase Voltage 13 

5. Modulated Phase Voltage 14 

6. Phase Current - 4500 RPM, No Modulation 20 

7. Fourier Spectrum of Phase Current - 4500 RPM, No Modulation 21 

8. Neutral Current - 4500 RPM, No Modulation 22 

9. Developed Torque - 4500 RPM, No Modulation 23 

10. Fourier Spectrum of Developed Torque - 4500 RPM, No Modulation 24 

11. Phase Current - 4500 RPM, Fifth and Seventh Harmonic Eliminated 26 

12. Fourier Spectrum of Phase Current - 4500 RPM, Fifth and Seventh 

Harmonics Eliminated 27 

13. Neutral Current - 4500 RPM, Fifth and Seventh Harmonics 

Eliminated 28 

14. Developed Torque - 4500 RPM, Sixth Harmonic Eliminated 29 

15. Fourier Spectrum of Developed Torque - 4500 RPM, Sixth Harmonic 

Eliminated 30 

16. Phase Current - 4500 RPM, Fifth, Seventh, Eleventh, and 

Thirteenth Harmonics Eliminated 32 

17. Fourier Spectrum of Phase Current - 4500 RPM, Fifth, Seventh, 

Eleventh, and Thirteenth Harmonics Eliminated 33 

18. Neutral Current - 4500 RPM, Fifth, Seventh, Eleventh, and 

Thirteenth Harmonics Eliminated 34 

19. Developed Torque - 4500 RPM, Sixth and Twelfth Harmonics 

Eliminated 35 

20. Fourier Spectrum of Developed Torque - 4500 RPM, Sixth and 

Twelfth Harmonics Eliminated 36 

21. Phase Current - 2250 RPM, No Modulation 37 

22. Fourier Spectrum of Phase Current - 2250 RPM, No Modulation 38 

23. Developed Torque - 2250 RPM, No Modulation 39 

24. Fourier Spectrum of Developed Torque - 2250 RPM, No Modulation 40 

25. Phase Current - 2250 RPM, Fifth and Seventh Harmonics Eliminated 42 

26. Fourier Spectrum of Phase Current - 2250 RPM, Fifth and Seventh 

Harmonics Eliminated 43 

27. Developed Torque - 2250 RPM, Sixth Harmonic Eliminated 44 

28. Fourier Spectrum of Developed Torque - 2250 RPM, Sixth Harmonic 

Eliminated 45 

29. Phase Current - 2250 RPM, Fifth, Seventh, Eleventh, and 

Thirteenth Harmonics Eliminated 46 

30. Fourier Spectrum of Phase Current - 2250 RPM, Fifth, Seventh, 

Eleventh, and Thirteenth Harmonics Eliminated 47 

31. Developed Torque - 2250 RPM, Sixth and Twelfth Harmonics 

Eliminated 48 

32. Fourier Spectrum of Developed Torque - 2250 RPM, Sixth and 

Twelfth Harmonics Eliminated 49 






ww w www iwwmv»wi am vmwinv wv. r? ■.» j ■ » ■ «t? ■ -% % ■ • ^ 


LIST OF ILLUSTRATIONS (CONTINUED) 

FIGURE PAGE 

33. Phase Current - 2250 RPM, a - 15°, Fifth, Seventh, Eleventh 

and Thirteenth Harmonics Eliminated 51 

34. Fourier Spectrum of Phase Current - 2250 RPM, a - 15°, Fifth, 

Seventh, Eleventh, and Thirteenth Harmonics Eliminated 52 

35. Developed Torque - 2250 RPM, a - 15°, Sixth and Twelfth 

Harmonics Eliminated 53 

36. Fourier Spectrum of Developed Torque - 2250 RPM, o » 15°, 

Sixth and Twelfth Harmonics Eliminated 54 

37. Phase Current - 2250 RPM, o - 40°, Fifth, Seventh, Eleventh, 

and Thirteenth Harmonics Eliminated 55 

38. Fourier Spectrum of Phase Current - 2250 RPM,a - 40°, Fifth, 

Seventh, Eleventh, and Thirteenth Harmonics Eliminated 56 

39. Developed Torque - 2250 RPM, a « 40°, Sixth and Twelfth 

Harmonics Eliminated 57 

40. Fourier Spectrum of Developed Torque - 2250 RPM, a - 40°, 

Sixth and Twelfth Harmonics Eliminated 58 

41. Phase Current - 450 RPM, No Modulation 60 

42. Fourier Spectrum of Phase Current - 450 RPM, No Modulation 61 

43. Developed Torque - 450 RPM, No Modulation 62 

44. Fourier Spectrum of Developed Torque - 450 RPM, No Modulation 63 

45. Phase Current - 450 RPM, Fifth and Seventh Harmonics Eliminated 65 

46. Fourier Spectrum of Phase Current - 450 RPM, Fifth and Seventh 

Harmonics Eliminated 66 

47. Developed Torque - 450 RPM, Sixth Harmonic Eliminated 67 

48. Fourier Spectrum of Developed Torque - 450 RPM, Sixth Harmonic 

Eliminated 68 

49. Phase Current - 450 RPM, Fifth, Seventh, Eleventh, and 

Thirteenth Harmonics Eliminated 69 

50. Fourier Spectrum of Phase Current - 450 RPM, Fifth, Seventh, 

Eleventh, and Thirteenth Harmonics Eliminated 70 

51. Developed Torque - 450 RPM, Sixth and Twelfth Harmonics 

Eliminated 71 

52. Fourier Spectrum of Developed Torque - 450 RPM, Sixth and 

Twelfth Harmonics Eliminated 72 

B.l Flow Chart for Selection of Candidate Modulation Functions 83 

B.2 Flow Chart for Optimizing Modulation Function 99 

B.3 Form of Modulation Functions 100 


vi 





LIST OF TABLES 


TABLE 


PAGE 


1 Motor Parameters 7 

2 Source Characteristics 7 


vii 










TrrTTgM ' J.T * y'\ ’> U' VlTl«IT. * '■ ■ ■■- ■ ■ ■ ■ A ' .T ■ ■ j. ’q 


SECTION I 
INTRODUCTION 


1.0 BACKGROUND 
1.1 Orientation 

Because of their inherent flexibility of control, reliability, and 
efficiency, uses of electromechanical actuators and electric drives are 
finding increased interests in aircraft applications. Three areas of 
technology advancement over the past decade are largely responsible for 
placing electromechanical energy converters into a favorable position 
when compared to hydraulic motors or actuators: 


(a). 

Development of rare earth permanent magnet (PM) motors with 
inherent high efficiency and high power to weight ratios. 

(b). 

Advancements in power level solid-state devices with high 
switching speeds. 

(c). 

Emergence of microprocessors allowing control capability 
and sophistication far surpassing that of analog systems 
while reducing the volume occupied. 

The conventional or commutator DC machine performance characteristics in 
the areas of speed control and position control are highly desirable. 

The brushless DC machine has vide range torque-speed characteristics like 
unto the commutator DC machine without the commutator-brush maintenance 
problems. In addition, the brushless DC machine with a permanent magnet 
rotor has certain superior features to the conventional DC machine: 

- 

Field excitation is eliminated, which removes the complexity 
of supplying power to a rotating member. Also, machine 
efficiency is increased due to absence of field excitation 
losses. 

(b). 

Higher speed design is possible for PM rotors than is feasible 
with wound rotors permitting Increased gear ratios, which leads 
to substantial reduction in electric machine power to weight 
ratios. 

(c). 

Thermal transfer characteristics are improved since the bulk 
of the losses (ohmic and core losses) are generated within the 
stationary member, allowing efficient implementation of fluid 
cooling. 
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1.2 Identification of Problem 

Brushless DC motor performance is reported in the literature, but it 
concentrates on the nature of instantaneous voltage and current along with 
average values of developed torque tl-93. These works typically describe 
systems that do not operate continuously at near zero speed, and thus, 
only average value of torque is of concern. However, the brushless DC 
machine inherently has oscillatory components of instantaneous torque 
at all frequencies that are integer multiples of six times the electrical 
radian of the stator impressed voltage. Since stator frequency is directly 
related to rotor position, these oscillatory components are in the fre¬ 
quency range of mechanical system response at low speed values. 

Published works relating to performance of brushless DC machines 
have not analyzed cases of position or sustained near zero speed operation; 
thus, pulsating torque components have been considered of no significant 
consequence. However, Williamson et al CIO] have mentioned that torque 
capability deteriorates at low speeds, and Demerdash and Nehl have shown 
some Instantaneous developed torque and power wave forms Cll3 without 
comment on the oscillatory components. 

Widespread usages of brushless DC motor drive systems in low-speed 
and position-control actuator applications are contingent upon development 
of control philosophies and hardware realizations of power conditioning 
arrangements that allow bidirectional power flow while resulting in 
Instantaneous motor developed torques that are free of harmonics in the 
range of response for coupled mechanical loads. 

1.3 Consideration of Alternative Technologies 

There are two basic power electronic circuits that are used as power 
conditioning links to couple brushless dc motor drives to a high frequency, 
three-phase AC source: 

(a) DC link inverter 

(b) Cycloconverter link 

Power circuits of these two basic approaches are illustrated in functional 
block form by Figure 1. Variations of each arrangement are made .according 
to the control philosophy adopted to satisfy required motor performance. 
Further, for low voltage (500 V or less), low current (200 A or less) 
applications, power circuits can be synthesized with transistors for 
controlled switching elements, while for high power level applications 
silicon controlled rectifiers (SCRs) must be used for controlled switching 
elements. For this study, SCR switching elements are assumed. 

1.3.1 Operation for Speeds Not Hear Zero . For the case of a DC link 
inverter driving a brushless DC motor that does not operate at near zero 
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Figure 1. Alternative Power Circuits 

(a) DC Link Inverter 

(b) Cycloconverter Link 
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speed, the rectifier may be a diode bridge while the chopper is used to 
vary magnitude of DC voltage applied to the inverter input terminals; or, 
the chopper may be eliminated and a phase-controlled converter used as 
the rectifier to vary inverter input voltage. 

For the case of a cycloconverter, the necessity of a DC link does 
not exist, and thus, there is only one power conditioning module inter¬ 
connecting the three-phase AC source with the brushless DC motor. 

Generally, the cycloconverter may be realized as either a midpoint or a 
full-bridge arrangement and either phase-control or synchronous envelope 
operation implemented. 

1.3.2 Operation for Speeds Near Zero . For near zero speed operation of 
a DC link inverter for high power level systems, the counter emf of each 
motor phase is small enough so that natural commutation of the inverter 
SCRs is prohibited. In such a case, the chopper is used to reduce 
inverter terminal voltage to zero sufficiently long to accomplish commu¬ 
tation of the inverter SCRs C3], Since the chopper must be controlled 
so as to enhance inverter commutation, the rectifier must be a phase- 
controlled converter to permit necessary inverter input voltage magnitude 
control. 

For near zero speed operation of a phase-controlled cycloconverter 
using circulating current free mode, discontinuous load current tends to 
occur, leading to increase in load current harmonics. Continuous 
current can be restored by changing to circulating current mode of control 
at the expense of increased losses due to circulating reactive current Cl2], 
Control is generally simplified at fractional hertz frequency operation by 
use of synchronous envelope control with circulating current free mode. 

1.3.3 Operation for Regenerative Power Flow . For the case of a DC link 
inverter, regenerative power flow requires further power circuit modifica¬ 
tion. Gating signals to inverter SCRs are suppressed and the associated 
shunting diodes function to form a full-bridge diode rectifier. A 
switching circuit is introduced in the DC link to reverse polarity of the 
voltage appearing at the DC terminals of the rectifier module. The 
rectifier must be a phase-controlled converter operated in the synchronous 
inversion mode. Reversal of motor developed torque requires three 
coordinated control actions: phase forward and suppression of inverter 
SCRs; polarity reversal switch activation; and, phase forward of phase- 
controlled converter SCRs. 

Regenerative power flow using the cycloconverter link is Introduced 
by simply delaying SCR firing angles beyond 90°. Rate of change of motor 
developed torque from a positive to a negative value is controlled by the 
rate at which the SCR firing angles are changed. Thus, the nature of 
transition from motoring to regeneration is determined by a single control 
action leading to smooth change with minimum time delay. 


4 








2.0 OBJECTIVE 


Of the presently available power conditioning technologies for 
linking the brushless DC motor to a high frequency polyphase AC source in 
a near zero speed or actuator application, the cycloconverter using syn¬ 
chronous envelope control with circulating free current mode of operation 
appears most promising in that it requires no external commutation cir¬ 
cuitry and it can be smoothly and quickly changed from motoring to 
regeneration by a single control action. This arrangement using a mid¬ 
point, three-pulse cycloconverter is chosen for study with pulse width 
modulation (PWM) control to eliminate those harmonics from the output 
waveform that lead to harmonics in the developed torque in the range to 
which coupled mechanical loads can respond. 

3.0 SCOPE 

The following tasks were carried out to accomplish the objective: 

(a) . Determine an appropriate mathematical model of the system for 

near zero speed operation. 

(b) . Determine the current harmonics to be eliminated in order to 

remove undesireable pulsations in motor developed torque. 

(c) . Develop a modulation function to use in control of phase 

voltage for eliminating the undesired current harmonics. 

(d) . Calculate motor performance over a wide speed range with and 

without the PWM phase voltage. 

(e) . Compare ohmic losses with and without PWM to assess the effect 

of harmonic elimination on efficiency. 

















SECTION II 


PROCEDURE 


1.0 SYSTEM DESCRIPTION 

Typical values for motor parameters were obtained by ratio of known 
values reported in the literature Cll] where the results are given in 
Table 1. The motor is wye connected with a maximum design speed of 
45000 rpm. 


TABLE 1. MOTOR PARAMETERS 
Parameter Symbol Value 


No. poles 
EMF constant 
Stator resistance 
Stator leakage Inductance 


P * 

K 0.0225 V-s/rad 

R_ 0.4 £2 

L a 25 yH 


Characteristics of the balanced three-phase source are listed in 
Table 2. 


TABLE 2. SOURCE CHARACTERISTICS 


Quantity 


Value 


Phases 3 

Line voltage 138 V 

Frequency 7950 Hz 


The power circuit schematic of a midpoint, three-pulse cycloconverter 
linking the high frequency, three-phase AC source to the brushless DC 
motor is depicted by Figure 2. Use of the switch in the midpoint or 
neutral line will be discussed later. 

2.0 ASSUMPTIONS 

When analyzing rare earth permanent magnet machines with non-magnetic 
retaining rings, it has been found that position dependence and interphase 
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coupling of machine inductances can be neglected Eli]. Such an approxima¬ 
tion leads to decoupled equations that can be used in networks formed by 
addition of power conditioning circuitry with minimum difficulty. Motor 
counter emfs are taken to be sinusoidal. Further, since the study of 
this report is concerned with low speed operation, the commutation overlap 
is very small compared with a period of the motor current wave form and 
can be neglected. 

The SCRs are modelled by a 0.02 ohm forward resistance (2 V at 100 A) 
and a blocking resistance described by 1000 i 2 ohms. 

3.0 MATHEMATICAL MODEL 

Referring to Figure 2, using the above assumptions, and applying 
Kirchhoff’s voltage law results in the following set of simultaneous 
differential equations to describe the motor electrical performance: 


di3_ 

dt 


r- <-M + T cn - e 3> 


(3) 


Ri, R 2 » and R 3 are the sum of the motor phase resistance and the resistance 
of the particular SCR between the source and phase of the motor at the 
particular Instant of time under analysis. The form of phase voltages 
< v an* v bn» v cn^ wil1 be dl 8 C“ 8 sed later. 

Neutral connection current is given by Kirchhoff’s current law as 


inN - ii + i 2 + h <*> 

As long as motor speed (u> m ) remains non zero, the electromechanical 
developed torque is given by 


T d - (ejij + e 2 i 2 + e 3 i 3 )/ Um (5) 

Equations (1) - (5) completely describe the instantaneous performance 
of the brushless DC motor. 










■T’.i " vw t <. ii i, t' / 






4.0 CURRENT HARMONICS TO BE ELIMINATED 

For both the case of a DC link inverter drive and the case of a 
synchronous envelope cycloconverter drive, the motor phase currents 
approach the pulse width controlled square wave of Figure 3. A Fourier 
series representation of the wave form shows that all odd, nontriplen 
harmonics exist in the phase current. 

i. ■ — 7 — cos (—) sin nut, n » 1,5,7,11,13, •• • (6) 

1 n n 6 



Figure 3. Approximate Form of Motor Phase Current 

Similarly, the remaining balanced phase currents are given by 

1 2 - *1 l I cos (Hp sin n(a.t - JiL.) (7) 

1 3 - I ^ cos (HI) sin n(ut + -^L.) (8) 

n b 3 

n - 1,5,7,11,13,••• 

Since the motor counter emfs (e 1 ,e 2 ,e.) are sinusoids of fundamental 
frequency and form a balanced three-phase set, each can be expressed as 
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E m sin («t + $) 


(9) 


e 2 - \ sin (ut + 4 - (10) 

e 3 “ Ejj, sin (wt + ♦ + -^-) (11) 

If equations (6) - (11) are substituted into (5), the simplified 
result for instantaneous developed motor torque is found. 

T d - £ Em 1 l h COB <T") {cosC(n-l)J u,t - <fr] - cosC(n+l)f ut + ♦]}(12) 
ir Q “ o o 6 

n « 1,5,7,11,13,••• 


From (12), it is seen that the instantaneous developed motor torque 
is made up of a constant term plus all multiples of the sixth harmonic 
of motor phase current frequency. Further, the multiples of sixth 
harmonic components of torque are a direct result of the current harmonics 
that lie immediately on either side of multiples of six. Specifically, 
the sixth harmonic of torque results from the fifth and seventh harmonic 
of current; and, the twelfth harmonic of torque results from the eleventh 
and thirteenth harmonic of current; etc. Clearly, the motor developed 
torque will be constant if the harmonics of motor phase current described 
by 6n ± 1, n - 1,2,3,••• are eliminated. 

The above discussed harmonics only need to be eliminated if their 
existence results in torque harmonics in the range of response of a 
coupled mechanical load. Since typical mechanical loads exhibit negligible 
response to frequencies above 20 hertz, it should only be necessary to 
eliminate those nontriplen, odd current harmonics beyond the first for 
which 


M < 


80 x 
6 P“ m 


(13) 


where p - number of poles and u m is the motor speed in rad/s. 

5.0 MODULATION OF PHASE VOLTAGE 

The cycloconverter system is voltage excited; thus, if the (n ± 1) 
harmonics are to be eliminated from motor phase currents, they must be 
eliminated from the phase voltage set. 

5.1 Development of Modulation Function 

Modulation function h(ut) that contains no selected harmonics up 
through M is derived in Appendix B. The balanced three-phase voltage set 
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to be applied for equations (1) - (3) is given by 


v an 

- v d h(ut) 

(14) 

v bn 

2 it 

- v d h(ut - — ) 

(15) 

v cn 

- v d h(u)t + ) 

(16) 


If is the radian frequency of the three-phase AC source, then v d is a 
periodic wave form described by 


v d (t) - 


and, v d (t + T) ■ 


where T 


It is shown in Appendix A that if the modulated phase voltages 
(van»vbn» v cn) are to contain none of the selected harmonics through M, 
it is necessary that the frequency of v<j be much greater than the 
frequency of h(t); however, this condition is easily satisfied at near 
zero speed. 

Figure 4 illustrates the unmodulated phase voltage. The modulated 
phase voltage is shown by Figure 5. 

5.2 Minimization of Slack Variable 

Development of the modulation function for use with the synchronous 
envelope cycloconverter in Appendix A required an extension of the well- 
known pulse width modulation techniques for elimination of harmonics in 
inverter drives Cl4] to include a slack variable. The modulation function 
contained no selected harmonics through M as determined by (13). For 
the reported work, the slack variable is always chosen so as to minimize 
the sum of the squares of the amplitude of the next higher (6n ± 1) 
current harmonic pair beyond harmonics selectively eliminated. Listings 
of the FORTRAN programs used to compute the modulation function are 
presented in Appendix B. 


V m sinCo^t + g- + o g ) , 0 < u^t < — (17) 

v d (t) (18) 

2ir 

3«L 
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5.3 Cycloconverter Connection Method 


Since the modulation function h(wt) does not repeat itself every ir/6 
radians on the interval Cir/6 , 5it/ 6D, attention must be directed to the 
cycloconverter connection scheme utilized so that the eliminated harmonics 
are not inadvertently reintroduced. Specifically, a conduction path must 
be provided so that all harmonics remaining in the phase voltage can flow 
if the phase current is to be of the same harmonic content. Two connection 
realizations are possible wherein no restriction on phase current harmonics 
are imposed by Kirchhoff’s current law. The first of these acceptable 
connections is that of the midpoint cycloconverter with the neutral lead 
in place. A second realization would be the full-bridge cycloconverter 
with total phase isolation. Since this latter connection requires 36 SCRs, 
it was not used in the reported study due to a factor of three on component 
requirement over the midpoint cycloconverter case. However, it should be 
pointed out that the full-bridge cycloconverter with its six-pulse output 
has a smaller ripple magnitude and may offer a slight efficiency advantage. 

5.4 Fourier Spectrum Analysis 

Nature of the waveforms for current and instantaneous developed 
torque are such that visual inspection to determine harmonic content is 
not practical or accurate. The listing of a FORTRAN program to find the 
Fourier coefficients is given in Appendix D. 

This program was used to generate the data for the normalized Fourier 
spectrum plots of Section III for harmonics up through forty-eight in order 
to verify thap attempt to selectively eliminate harmonics was successful. 




















SECTION III 
RESULTS 


1.0 SUMMARY 

Three speed conditions were selected for numerical study to evaluate 
the effectiveness of the harmonic elimination technique in removing 
undesired torque pulsations from the midpoint cycloconverter driven 
brushless DC motor system when using synchronous envelope control: 

(a) . 4500 rpm (10% speed case) 

(b) . 2250 rpm ( 5% speed case) 

(c) . 450 rpm ( 1% speed case) 

Although these speeds are great enough that elimination of harmonics as 
determined by equation (13) would not be necessary, there are two reasons 
for their selection. The first is conservation of computer time. Computa¬ 
tion of a set of data for the 450 rpm (1% speed case) required approxi¬ 

mately 8 hours of CPU time on a Control Data PDP 11/44 computer. Since 
many trial-and-error runs were necessary for each data set, runs at lower 
motor speed become prohibitive timewise. Secondly, if the harmonics can 
be successfully eliminated at higher speeds, then there is not difficulty 
to be encountered at reduced speed. Justification lies in the fact that 
the undeslred harmonics are totally eliminated from the phase voltage set. 
Any re-introduction of the undesired harmonics into the phase current 
wave forms would be attributable to commutation overlap which is dependent 
on the L a /R time constant that determines current decay at the end of 
each phase current conduction sequence. Since this time constant is not 
speed dependent, the commutation overlap interval is a larger portion of 
the current wave form period at high speed than it is at low speed. Thus, 
the departure of phase current from the form of the phase voltage is 
greater at high speed than at low speed. It is concluded that if 
undesireable harmonics are not re-introduced at high speed, then success 
at low speed is guaranteed. 

For each speed case, the nature of current and torque were examined 
for application of a balanced set of phase voltages of the following 
nature: 

(a) . No harmonic eliminated. 

(b) . Fifth and seventh harmonics eliminated. 

(c) . Fifth, seventh, eleventh, and thirteenth harmonics eliminated. 
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A brushless DC machine system has an extra Input or degree of freedom 
over either a conventional DC machine or a synchronous machine drive 
system in that the angular displacement between the stator and rotor 
magnetic fields can be specifically controlled. For this study, this 
extra degree of control freedom is implemented by specification of the 
angle (a) measured from the onset of positive phase current conduction 
to the positive going zero crossing of the associated motor phase counter 
emf. A convenient value (a - 30°) is used for most of the study; however, 
since selection of a shifts the position of motor phase counter emf with 
respect to phase current, the nature of conamitation overlap can be altered 
by the value of a. The effect of changing a over the range from 15° to 
40° is studied for one speed case. 

2.0 PERFORMANCE STUDIES 


Equations (1) - (3) can be written in compact matrix form as 


dfi “ ^i+^-JLe 


( 20 ) 



'ir 

(21) ; e - laK 
“ P 

sin(ujt-oX 

1 - 

i 2 

sin(ut-a- IlL ) 

3 


_ i 3 . 


_sin(uit-a+ ^ ) 


( 22 ) 


CR] - 

" R 1 

r 2 

(23) ; 24 - v d (t) 

h(wt) 

h(u»t - ^ ) 


R 

L 3 J 


h(ut + ^ ) 


(24) 


Equation (20) is a nonlinear set of differential equations due to the 
elements of CR] being functions of current. Any attempt at linearization 
for a wide range of i_would unacceptably alter the physical nature of 
the problem by removing the rectifying characteristic of the SCRs. Conse¬ 
quently, a closed form solution of (20) is not possible and a numerical 
integration must be implemented. A fixed increment (1.0 x 10~*>s) fourth- 
order Runge-Kutta numerical integration method is utilized for this work. 

After forming a modulation function through use of the programs of 
Appendix B, selecting a particular motor speed (u^ - 2u/p), and arbitrarily 
choosing a value of SCR firing angle a g (see equation (17)), numerical 
solution to (20) is implemented subject to initial conditions i_(0) ■ 0. 
Integration is continued until steady-state is reached. Instantaneous 
values of developed torque are calculated according to (5) and numerically 









averaged over the last steady-state cycle of integration to give the 
average value of developed torque. This value of average torque is 
compared with a specified value for the operating point. If the average 
value calculated is different from the specified value, then ot g is 
appropriately adjusted and the numerical integration is repeated until the 
difference magnitude between calculated and specified torque are suffi¬ 
ciently close in an epsilon sense. For all cases studied, the specified 
torque was taken as 1.5 N-m. 

The listing of a FORTRAN program to compute the performance as 
described above is presented in Appendix D. In addition to the discussed 
computations, the program also calculates the RMS value of phase current 

over the last steady-state cycle as ( i, £ i 2 At)**. Also, values of phase 

current and instantaneous torque are stored over the last steady-state 
cycle for use by the Fourier Spectrum Program of Appendix C. 

Results of the various speed cases and harmonic elimination studies 
are summarized in the paragraphs that follow. 

2.1 4500 RFM (10% Speed Case) 

2.1.1 Unmodulated Phase Voltage . 

Input Data 

Motor speed, - 4500 rpm 

Motor frequency, f - 150 Hz 
SCR firing angle, a s - 58.68° 

Specified average torque, T^^ - 1.5 N-m 


Output Data 

Average torque, Tdav “ 1.501 N»m 
RMS phase current, I^g - 72.77 A 


A plot of phase current ij is given by Figure 6. The associated 
Fourier spectrum of ij is shown by Figure 7. Predominantly triplen 
harmonic current that flows through the neutral connection is depicted by 
Figure 8. 

A plot of the instantaneous developed torque where the sixth harmonic 
is obvious is given by Figure 9. A Fourier spectrum of the instantaneous 
torque is shown in Figure 10 where existence of all multiples of six as 
theoretically predicted are noted. 
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Figure 7. Fourier Spectrum of Phase Current - 4500 RPM, No Modulation 

















in 

n 


in 


r<; 




IO 




23 
















Figure 10. Fourier Spectrum of Developed Torque - 4500 RPM, No Modulation 







2.1.2 Elimination of Sixth Harmonic of Torque . The sixth harmonic of 
torque Is eliminated by removing the fifth and seventh harmonics of phase 
current. 


Motor speed, - 4500 rpm 

Motor frequency, f » 150 Hz 
SCR firing angle, o g - 46.45° 

Specified average torque, T dav - 1.5 N*m 
Angles of modulation function, o. ■ 40.76° 
ai - 47.73° 
o 3 - 58.65° 


Average torque, T, - 1.5001 N*m 

RMS phase current,^KMS " 7 **24 A 


A plot of phase current 1^ is given by Figure 11. The associated 
Fourier spectrum of i^ is shown by Figure 12, where negligibly small fifth 
and seventh harmonics are noted. Current flow through the neutral line 
is depicted by Figure 13. 

A plot of the instantaneous developed torque is given by Figure 14. 
The Fourier spectrum of the developed torque is shown in Figure 15 where 
it is seen that the sixth harmonic is near zero. 

2.1.3 Elimination of Sixth and Twelfth Harmonics of Torque . Both the 
sixth and twelfth harmonics of torque are eliminated by removing the fifth, 
seventh, eleventh, and thirteenth harmonics of phase current. 


Motor speed, * 4500 rpm 

Motor frequency, f - 150 Hz 
SCR firing angle, o # - 42.63° 

Specified average torque, T<j av * 1«5 N»m 
Angles of modulation function, ai ■ 38.73° 
a 2 - 42.13° 
o 3 - 57.25° 
ct 4 - 61.93° 
u 5 - 66.86° 


Average torque, T<j av - 1.505 N.m 

RMS phase current, - 74.59 A 
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Figure 12. Fourier Spectrum of Phase Current - 4500 RPM, Fifth and 
Seventh Harmonics Eliminated 
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4500 RPM, Fifth and Seventh Harmonics Eliminated 
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A plot of phase current ij^ Is given by Figure 16. The associated 
Fourier spectrum of ij is shown by Figure 17 where negligibly small fifth, 
seventh, eleventh, and thirteenth harmonics are noted. Current flow 
through the neutral line is depicted by Figure 18. 

A plot of the instantaneous developed torque is given by Figure 19. 
The Fourier spectrum of the developed torque is shown in Figure 20 where 
it is seen that the sixth and twelfth harmonics are near zero. 

2.2 2250 RPM <5% Speed Case) 

2.2.1 Unmodulated Phase Voltage . 

Input Data 

Motor speed, • 2250 rpm 

Motor frequency, f - 75 Hz 
SCR firing angle, 8 - 61.93° 

Specified average torque, T dav - 1.5 N*m 

Output Data 

Average torque, T dav - 1.499 N-m 

RMS phase current, IrmS “ 69-53 A 


A plot of phase current ^ is given by Figure 21. The associated 
Fourier spectrum of ij is shown by Figure 22. 

A plot of the instantaneous developed torque is given by Figure 23. 

A Fourier spectrum of the developed torque is shown by Figure 24. 

2.2.2 Elimination of Sixth Harmonic of Torque . The sixth harmonic of 
torque is eliminated by removing the fifth and seventh harmonics of phase 
current. 


Input Data 

Motor speed, n m ■ 2250 rpm 

Motor frequency, f - 75 Hz 
SCR firing angle, a s - 50.21° 

Specified average torque, T^av * 1-5 N*m 
Angles of modulation function, a. - 40.76° 
a 2 - 47.73° 
a 3 - 58.65° 


Output Data 

Average torque, T da ^ - 1.5001 N-m 
RMS phase current, Irms " 73.96 A 
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Figure 16. Phase Current - 4500 RPM, Fifth, Seventh, Eleventh 
and Thirteenth Harmonics Eliminated 













Figure 17. Fourier Spectrum of Phase Current - 4500 RPM, Fifth, 

Seventh, Eleventh, and Thirteenth Harmonics Eliminated 
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A plot of phase current 1^ is given by Figure 25. The associated 
Fourier spectrum is shown by Figure 26. 

A plot of the instantaneous developed torque is given by Figure 27. 

The Fourier spectrum of the developed torque is shown by Figure 28. 

2.2.3 Elimination of Sixth and Twelfth Harmonics Of Torque . Both the 
sixth and twelfth harmonics of torque are eliminated by removing the fifth, 
seventh, eleventh, and thirteenth harmonics of phase current. 

Input Data 

Motor speed, n,,, ■ 2250 rpm 

Motor frequency, f ■ 75 Hz 
SCR firing angle, a a * 47.31° 

Specified average torque, T dav - 1.5 N-m 
Angles of modulation function, oi * 38.73° 
a 2 - 42.13° 
a 3 - 52.25° 
a. - 61.93° 
a 5 - 66 . 86 ° 

Output Data 

Average torque, T dav - 1.5003 N-m 
RMS phase current, Irms * 74.64 A 

A plot of phase current i^ is given by Figure 29. The associated 
Fourier spectrum of i^ is shown by Figure 30. 

A plot of the instantaneous developed torque is given by Figure 31. 

The Fourier spectrum of the developed torque is shown by Figure 32. 

2.2.4 Elimination of Sixth and Twelfth Harmonics of Torque (a " 15°) . The 
angle from the onset of positive phase current to the positive going 
crossing of the associated motor phase counter emf is changed from 30° to 
15° to examine effect on current harmonics. 


Input Data 


Motor speed, tig, ■ 2250 rpm 

Motor frequency, f • 75 Hz 
SCR firing angle, a s - 59.84° 
Specified average torque, T<| av * 
Angles of modulation function, ai 
°2 
°3 
a 4 
°5 


1.5 N*m 

- 38.73° 

- 42.13° 

- 52.25° 

- 61.93° 

- 38.73° 
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Figure 26. Fourier Spectrum of Phase Current - 2250 RPM, Fifth 
and Seventh Harmonics Eliminated 
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Fourier Spectrum of Developed Torque - 2250 RPM, 
Sixth Harmonic Eliminated 
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Fourier Spectrum of Developed Torque - 2250 RPM, 
Sixth and Twelfth Harmonics Eliminated 





Average torque, T dav * 1.502 N*m 
RMS phase current, “ 53.33 A 


A plot of phase current i, is given by Figure 33. The associated 
Fourier spectrum Is shown by Figure 34. 

A plot of Instantaneous developed torque Is given by Figure 35. 

The Fourier spectrum of the developed torque is shown by Figure 36. 

The re-lntroduced harmonic magnitudes are reduced to a lower level 
than for a ■ 30°. 

2.2.5 Elimination of Sixth and Twelfth Harmonics of Torque (a ■ 40°) . 
The angle from onset of positive phase current to the positive going 
crossing of the associated motor phase counter emf is set at 40° to 
examine the effect on current harmonics. 

Input Data 

Motor speed, n m - 2250 rpm 

Motor frequency, f « 75 Hz 
SCR firing angle, u 8 - 16.77° 

Specified average torque, T dav - 1.5 N-m 
Angles of modulation function, 04 - 38.73° 
o 2 - 42.13° 

03 - 52.25° 
o 4 - 61.93° 

05 - 66 . 86 ° 


Output Data 


Average torque, T dav - 1.5003 N«m 
RMS phase current, IrmS " 107.52 


A plot of phase current i 1 is given- by Figure 37. The associated 
Fourier spectrum of phase current is shown by Figure 38. 

A plot of instantaneous developed torque is given by Figure 39. 

The Fourier spectrum of developed torque is shown by Figure 40. 

Although the angle power factor angle (o + 30°) has been extended to 
75° resulting in a large current to maintain the same power flow, no 
significant re-introduction of undesired harmonics is noted. 
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Phase Current - 2250 RPM, a - 15°, Fifth, Seventh, 
Eleventh, and Thirteenth Harmonica Eliminated 








00 T 








































Degrees 

Figure 37. Phase Current - 2250 RPM, a - 40°, Fifth, Seventh, 
Eleventh, and Thirteenth Harmonics Eliminated 
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2.3 450 RPM (1% Speed Case) 

2.3.1 Unmodulated Phase Voltage . 

Input Data 

Motor speed, n m ■ 450 rpm 

Motor frequency, f ■ 15 Hz 
SCR firing angle, a g - 64.52° 

Specified average torque, T(j av * 1.5 N*m 

Output Data 

Average torque, T d ■ 1.4996 N*m 
RMS phase current, Ig^s * 67.1 A 


A plot of phase current ii is given by Figure 41. The associated 
Fourier spectrum is shown by Figure 42. 

A plot of Instantaneous developed torque is given by Figure 43. A 
Fourier spectrum of developed torque is shown by Figure 44. 

The .instantaneous wave form plots for the low speed case do not dis¬ 
play all of the rippled due to the lack of capacity of the plot routine 
to store sufficient points. 

2.3.2 Elimination of Sixth Harmonic of Torque . The sixth harmonic of 
torque is eliminated by removing the fifth and seventh harmonics of phase 
current. 


Input Data 

Motor speed, ■ 450 rpm 

Motor frequency, f ■ 15 Hz 
SCR firing angle, a g ■ 52.53° 

Specified average torque, T dflV - 1.5 N*m 
Angles of modulation function, a, ■ 40.76° 
a 2 - 47.73° 
a 3 - 58.65° 

Output Data 

Average torque, T<j av - 1.4996 N*m 
RMS phase current, Ipjjs ■ 75.2 A 


A plot of phase current i 1 is given by Figure 45. The associated 
Fourier spectrum is shown by Figure 46. 
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A plot of instantaneous developed torque is given by Figure 47. 
The Fourier spectrum of developed torque is shown by Figure 48. 

2.3.3 Elimination of Sixth and Twelfth Harmonics of Torque . Both the 
sixth and twelfth harmonics of torque are eliminated by removing the 
fifth* seventh* eleventh* and thirteenth harmonics of phase current. 

Input Data 


Motor speed, - 450 rpm 

Motor frequency, f ■ 15 Hz 
SCR firing angle, a g - 49.78° 
Specified average torque, T dav 
Angles of modulation function. 


- 1.5 N*m 


°1 

a 2 

°3 

“4 

®5 


38.73° 

42.13° 

52.25° 

61.93° 

66 . 86 ° 


Output Data 

Average torque, T dav - 1.507 N»m 
RMS phase current Irms - 77.38 A 


A plot of phase current i^ is given by Figure 49. The associated 
Fourier spectrum of phase current is shown by Figure 50. 




A plot of instantaneous developed torque is given by Figure 51. The 
Fourier spectrum of developed torque is shown by Figure 52. 
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Figure 46, 

















































10Q 



Harnonic Nanber 

Figure 50. Fourier Spectrum of Phase Current - 450 RPM, Fifth, 

Seventh, Eleventh, and Thirteenth Harmonics Eliminated 
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Developed Torque - 450 RPM, Sixth 
and Twelfth Harmonics Eliminated 
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SECTION IV 


DISCUSSION 


1.0 INCREASE IN OHMIC LOSSES 

Since modulation of the wave form definitely redistributes the 
harmonic content of the phase current wave form, an examination of the 
change in winding ohmic losses is in order. For the 4500 rpm case it 
is noted that the RMS phase current is found to be as follows: 

(a) . Irms " 72.77 A for no modulation 

(b) . Igjjg - 74.24 A for fifth and seventh harmonic eliminated 

(c) . IrmS " 74.59 A for fifth, seventh, eleventh, and thirteenth 

harmonic eliminated 

The increase in ohmic losses assuming a constant winding resistance is 
given by the ratio of RMS current for a changed condition to the RMS 
value of current for the unmodulated wave form. It is concluded that 
the ohmic losses are increased by 2.1% for the case of elimination of 
fifth and seventh harmonics. The ohmic losses are increased by 2.5% for 
the case of elimination of the fifth, seventh, eleventh, and thirteenth 
harmonics. 

Calculated data for the 2250 rpm and 450 rpm cases show increases 
of 7 to 15 Z; however, there is no theoretical basis for greater increase 
in these cases. The reason for this apparent larger increase for the 
lower speed cases is a numerical problem due to storing and calculating 
the RMS value based on a limited number of data points (2000) for trans¬ 
mittal to the plot routine. As a result, the high frequency ripple 
currents are sampled at irregular intervals giving rise to inaccuracy in 
calculation of the RMS value. For the case of 4500 rpm, the data point 
limit allows storage of several points during each ripple cycle, thus, the 
resulting RMS values are considerably more accurate than for the other 
two speed cases. 

2.0 CYCLOCONVERTER MODE CHANGE 

As discussed and illustrated by Figure 2, a neutral connection is 
necessary to successfully eliminate selected phase current harmonics. 

Since harmonic elimination is only necessary at speeds near zero, and 
since an ohmic loss penalty is incurred as a result, it might be 
deslreable to not operate with PWM control over the full range of speed. 
The Harmonic Eliminator Switch in the neutral line of Figure 2 could be 
opened for higher speeds and the cycloconverter mode of operation changed 
to phase-control from synchronous envelope. 
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3.0 EXTENSION TO HIGHER HARMONIC ELIMINATION 

For specific illustration, this study has only considered elimination 
of torque harmonics through the twelfth. As motor speed approaches zero, 
it may be necessary to eliminate even higher torque harmonics. Theoreti¬ 
cally, there is no limitation to be encountered. 
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SECTION V 


CONCLUSIONS AND RECOMMENDATIONS 


The work performed by this study shows that systematic reduction of 
harmonic torque pulsations in brushless DC drives is feasible by use of 
selective current harmonic elimination. Further, there appears to be 
only a small ohmic loss penalty. 

Further theoretical study should be made to properly assess quanti¬ 
tatively the effect of increased motor leakage inductance and of interphase 
reactors on "reintroducing" the harmonics eliminated from the phase 
voltage due to extending the commutation overlap. 

Also, a hardware realization of a midpoint cycloconverter brushless 
DC drive using the control principles formulated in this study should be 
made to verify the results. 
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APPENDIX A 


HARMONICS OF PWM TIME FUNCTION 


A modulation function h(wt) was derived in Appendix B and used in 
Section II - 5.0 to form a modulated phase voltage. For example. 

Van ■ h <“t) v d (t). 

Harmonic content of a modulated voltage is of concern. The study 
is best understood by use of a typical example. Let h(wt) contain all 
odd harmonics except the fifth and seventh, then 


h(u>t) * ajcos wt + a 3 cos 3ut + a g cos 9u>t + a 11 cos llwt + ••• 

+ SjjCos nut + ••• (A.l) 


Assume the v d (t) contains a DC term plus one high frequency harmonic, then 

v d (t) - b Q + b m cos mut (A.2) 

Forming v^ as the product of h (wt) and v d (t) and simplifying yields 

v d (t) - CjCOs ut + c 3 cos 3ut + c 9 cos 9ut + c 11 cos llwt + ••• 

+ d^jCosCm+lJwt + d^^coB (m-l)ut + d^cosOirt^wt 
+ d m _ 3 cos(m-3)u)t + ••• (A. 3) 


It is seen from (A.3) that through the cross product terms, the PWM 
expression may contain frequencies that were selectively eliminated from 
h(<i>t). However, the coefficients (d^i, d m _i, ...) of (A.3) are products 
of ajj and bn. Since an decreases as l/ n , if m is large, then the 
coefficients of the low frequency components of h(ut) resulting cross 
product terms will be negligibly small. It is concluded that the fre¬ 
quency of the AC source must be large compared with the motor frequency 
if the PWM function is to not contain the harmonics selectively eliminated 
from the modulating function h(ut). But, at near zero speed, the source 
frequency is large compared with the motor frequency. 
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MODULATION FUNCTION PROGRAMS 


1.0 INITIAL SOLUTION 

As first procedure, a candidate or initial modulation function, h(wt), 
is found to serve as a basis for an optimization search to determine the 
final modulation function. 

1.1 Theory of Harmonic Elimination 

Extending Patel and Hoft's pulse width modulation technique C14D to 
include a slack variable, it is possible to find a modulation function 
that has selected harmonics elimination while existing on the quarter wave 
interval from 30° to 90°. 

The Fourier series representation of the modulation function assuming 
odd quarter wave symmetry is 

f(ut) « £ ajjSinGiwt) , for odd n 

where 

a 2 a 4 a M 

ajj - / a sin(nut)dut + / sin(nut)dait + ••• sin (nut) dm t, (B.l) 


30° < a x < o 2 < < «„ < 90° 


This reduces to 


I 


*n " " k-1 ( " 1} COS<m k ) (B - 2) 

for M both odd and even. The number of harmonics eliminated is in general 
the same as the number of switching transitions per quarter cycle. For M 
greater than 2 the equations cannot be solved in closed form, so numerical 
methods are needed to obtain solutions. The following normalized system 
of equations is to be used to obtain solutions: 


f nl " I (-1)- cosGijafc) - 0 

f n 2 " I <-l) k+1 cos(n 2 a k ) - 0 


f nM “ I (-l) k+1 cos(n M a k ) » 0 










where are the harmonics to be eliminated. 


Many of the solutions of (B.3) are meaningless or impractical, so it 
is necessary to search by brute force using a loose convergence criterion 
for candidate solutions, and then to use a numerical method such as Newton 
Raphson to obtain a refined solution. 

The basic method of selecting candidate modulation function for 
further optimization is illustrated by the flow chart of Figure B.l. A 
candidate set is one that meets the following criteria: 


f squar ^ f min 

f squar " f nl + f n2 *** f nM 

If for the candidate modulation function, 30° < aj < a 2 < * * * < a M < 90°, 
then the solution is considered for optimization. Otherwise, the solution 
is discarded. 








SELECT M HARMONICS 
AND SLACK HARMONIC 
TO BE ELIMINATED 



Figure B.l. Flow Chart for Selection of Candidate 
Modulation Functions 












1.2 Initial Solution Program 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 
c main program c 


ccecceeccccccccecccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
on- number of harmonics to be eliminated 
c this includes the slack variable 

• this is the amount that each alpha 
will be incremented by during the 
infinite search 

- newton rapson method is used to 
search for a solution if the sum 
of squares is less than fmin 

- these are the harmonics to be reduced 

- this is the minimum value of alpha 
for the infinite search 

- not used 

- not used 

dimension a(10) t r(10),f(10) 
common pi f conv,iwr,del 
real inc 
call sopen(O.O) 
pi=3.141592653 


inc 


fmin 


r(i) 

amin 


iwr 

del 


klminsO 

k2min=0 

k3min=0 

k4min=0 

k5mins0 

k6min=0 

k1=1 


k2*1 

k3»1 

k4s1 

k5=1 

k6*1 

conva180.0/pi 

read(9,») n 

write(6 f ») 1 ns’.n 

read(9,*) Inc 

write(6,*) * inea'.inc 

read(9 f # ) fmin 

wrlte(6,*) 1 fmina',fmin 

read(9,») (r(i),lsi, 7 ) 

write(6,») * harmonicsa* t (r(i),ial,n) 

read(9, § ) amin 

write(6,*) ' amins'.amin 

read(9, # ) iwr 

wrlte(6,*) ' iwra»,iwr 











max=int( (90.0-amin)/inc ) 

k1max=int( (90.0-amin)/inc ) 

k2maxsmax 

k3max=max 

kUmaxsmax 

kSoaxsnax 

k6max=max 

isepsintC 1.0/lnc ) +1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

e 

c the following are essentially do loops but were 

c necessary to get past a problem with the compiler 

c 

klainsl 

100 klsklmin-l 

105 k1=k1+1 

if( n .eq. 1 ) go to 1000 

k2mln=k1+lsep 
200 k2=k2mln-1 

205 k2=k2+1 

lf( n .eq. 2 ) go to 1000 

k3min=k2+isep 
300 k3=k3mln-1 

305 k3*k3+1 

lf( n .eq. 3 ) go to 1000 

k4min=k3+isep 
400 k4=k4min-1 

405 k4sk4+1 

if( n .eq. 4 ) go to 1000 

k5mln*k4+lsep 
500 k5*k5min-1 

505 k5*k5+1 

lf( n .eq. 5 ) go to 1000 

k6ainxk5+lsep 
600 k6*k6mln-1 

606 k6*k6«-1 

IfC n .eq. 6 ) go to 1000 

1000 a(1M float(kl )*lnc+amln )/conv 
a(2)*( float(k2) # inc+amln )/conv 
a(3)«( float(k3)*lnc-famin )/conv 
a(4)s( float(k4)*inc-famln )/conv 
a(5)«( float(k5)*inc+amin )/conv 
a(6)«( float(k6)»inc+amln )/conv 
call ff( a,r,f,n ) 
oall square( f t fsquar,n ) 

lf( fsquar .It. fmln ) call look( a,r,f,fsquar,n ) 
if( k6 .It. k6max .and. n .ge. 6 ) go to 606 

lf( k5 .It. k5max .and. n .ge. 5 ) go to 505 

if( k4 .It. k4max .and. n .ge. 4 ) go to 405 

lf( k3 .It. k3max .and. n .ge. 3 ) go to 305 
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3000 

3100 

3200 


if( k2 .It. k2max .and. n .ge. 2 ) go to 205 
lf( kl .It. klmax .and. n .ge. 1 ) go to 105 
wrlte(6,3100) 

write(6,3000) ( a(i)*conv,i=1,n ) 

call sclose(O.O) 

format( 10<e15.7,5x> ) 

format(* a*') 

formate harmonicas') 

stop 

end 
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occccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cocccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine a c 

c c 

ccccccccceccececccccoccocccccccccccceccccccccccccccccccccccccccccccecccc 
cccccccccoceococcccecccccccccccccccecccccccccccccccccccccccccccccccccecc 
subroutine fa( a,da,n ) 
dimension a(10),da(10) 
common pi»conv,iwr,del 

do 200 lsl,n 
a(i)=a(i)-da(i) 

200 continue 
return 
end 








cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 
c subroutine close c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sclose(dummy) 
closed) 
close(2) 

close(3) 

closed) 
close(7) 
closed) 
closed) 
return 













cccccccccccccccccccccccccccccccccccccccccccccocccccccccccccccccccccccccc 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c c 

c subroutine f c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccocccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecccc 
subroutine ff( a,r,f,n ) 
dimension a(10),r(10),f(10) 
common pi,conv,iwr,del 
do 200 isl.n 
sign=1.0 
f(i)=0.0 

do 100 Js1,n 

f(i)sf(i)+sign*cos( r(i)«a(j) ) 
signs-sign 

100 continue 

200 continue 
return 
end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

c subroutine look c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccecccccccccoccccccccccccccccccccccccccccoeccccccccccccccccccccccccccc 
subroutine look( a,r,f,fsquar,n ) 

dimension da(10),r(10),a(10) f f(10),df(10,10),dfi(10,10) 

common pi,conv,iwr f del 

write(6,3050) 

write(6,3050) 

write(6,3100) 

write(6,3000) ( a(i)*conv,i=1,n ) 
write(6,3150) 

write(6,3000) ( f(i),i=1,n ) 

write(6,3200) 

write(6,3000) fsquar 

kount=30 

fminsO.01 

delta=0.001 


do 2500 loop=l.kount 
call fdf( a,r,f,df,n ) 
call matv( df,n,det,dfi ) 
call multK f,dfi,da,n,n ) 
call fa( a,da,n ) 
call ff< a,r,f,n ) 
call square( f,fsquar,n ) 

if( loop .gt. 5 .and. fsquar .It. fmin ) go to 2900 
2500 continue 

2900 write(6,3050) 

write(6,3100) 

write(6,3000) ( a(i)*conv,is1,n ) 
write(6,3150) 

write(6,3000) ( f(i),is1,n ) 
write(6,3200) 
wrlte(6,3OOO0 fsquar 
3000 formatC 10(e10.4,5x) ) 

3050 formate ' • ) 

3100 formate ’ a=' ) 

3150 formate ' f=' ) 

3200 formate ' fdquars’ ) 

4000 formate'test',13) 

5000 return 

end 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C SUBROUTINE MATV C 

c c 

c c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE MATV(A,N,DET,B) 

DIMENSION INDEX(20,2),IPVOT(20),A(10,10),B(10,10),PIVOT(20) 

ONES1.0 
ZEROsO.O 

c 

DO 10 1=1,N 

DO 5 J=1,N 
B(I,J)=A(I,J) 

5 CONTINUE 

10 CONTINUE 
DET=ONE 
DO 15 J=1,N 
15 IPVOT(J)=0 

DO 75 1=1,N 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 15 STATEMENTS C 
C FOR SEARCH FOR PIVOT ELEMENTS C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
TsZERO 
DO 40 Jsl.N 

IF(IPV0T(J)-1) 20,40,20 
20 DO 35 Ksl.N 

IF(IPV0T(K)-1) 25,35,100 
25 CONTINUE 

DlsABS(T) 

D2=ABS( B(J,K) ) 

IF(D1-D2)30,35,35 
30 IROWsJ 

ICOLsK 
TsB(J,K) 

35 CONTINUE 

40 CONTINUE 
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mm wans ul to w s n '.«■ k lv.v.w.*.".* 


IPVOT(ICOL)=IPVOT(ICOL)+1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 10 STATEMENTS PUT C 
C PIVOT ELEMENT ON DIAGONAL C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF(IROW-ICOL) 45,55,45 
45 DETs-DET 

DO 50 L=1,N 
TsB(IROW.L) 

B(IR0W,L)=B(IC0L,L) 

50 B(IC0L,L)=T 

55 INDEX(I,1)=IR0W 

INDEX(I,2)=IC0L 
PIV0T(I)=B(IC0L,ICOL) 

DET=DET*PIV0T(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 3 STATEMENTS TO DIVIDE PIVOT C 

C ROW BY PIVOT ELEMENT C 

C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

B (ICOL,ICOL)=0NE 
DO 60 L=1,N 

60 B(IC0L,L)=B(IC0L,L)/PIV0T(I) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"JCCCCCCCCCCCCCCCCCCCCCCCC 


c c 

C FOLLOWING 7 STATEMENTS TO REDUCE C 
C NON-PIVOT ROWS C 
C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 75 LIsl.N 
IF(LI-ICOL) 65,75,65 
65 T*B(LI,IC0L) 

B(LI,ICOL)=ZERO 
DO 70 L*1,N 

70 B(LI,L)aB(LI,L>-B(ICOL,L)»T 

75 CONTINUE 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

c FOLLOWING 11 STATEMENTS TO C 
C INTERCHANGE COLUMNS C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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DO 95 1=1,N 
LsH-I+1 

IF(INDEX(L,1)-INDEX(L,2) ) 85,95,85 
JROWsINDEX(L,1) 

JC0L=INDEX(L,2) 

DO 90 K*1,N 

T=B(K,JROW) 

B(K,JROW)=B(K,JCOL) 

B(K,JCOL)=T 
CONTINUE 
CONTINUE 
RETURN 
END 









cccccccceecccccccccccccccccccccccecccccccccccccccccccccccccccccccccccccc 

cocccccccccecccccccccccccccccccccccccccccccccccccccccccccccccceccccccccc 


c subroutine multiply c 

c c 

cccccccccececccccccececccocccccceccccccccccecccccccccccceccccccccccccccc 
cceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine multi( a t b,c,10,k0 ) 
dimension a(10),b(10,10),e(10) 
do 200 is1,10 
c(i)s0.0 
do 100 k=1,k0 

c(l)*c(i)+b(l,k)*a(k) 

100 continue 

200 continue 
return 
end 
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ccccccccccccccccccccccccceccccccccccceccccccccccccccccccccccccccccccoccc 
coeccccccccccccccccceccccecccccccccecceccccccccccccccccccccccccccceccccc 
c c 

c subroutine open c 

e c 

cceccccccccccccccecccccccccccccccccccccccccccccccccccccccccccccecccceccc 
eccccccceccccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sopen(dummy) 
open(unit=9,files »input") 
rewind(9) 
return 
end 














ceeccccccccccccececeeecccccccccccccccccccccccccccccccccccccccccccccccccc 
ceceeccccccccoccccccoccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine square c 

c c 

cccocceeeccccccceeeeeccccccccceccccccccocccccccccccccccccccecccccecccccc 
eccececeeeececeeoooooeeccocecocecccccccoeccecccccccccccccccccccccccccccc 
subroutine square( f,fsquar,n ) 
dimension f(10) 
fsquar=0.0 
do 200 i=1,n 

fsquar*fsquar♦f(i)*f(i) 

200 continue 
return 
end 
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2.0 OPTIMIZATION OF MODULATION FUNCTION 
2.1 Discussion of Procedure 

A candidate solution for the modulation function is selected from 
the output of the Initial Solution Program for optimization. The slack 
harmonic of the Initial solution is changed in small Increments and fed 
to a Newton-Raphson routine to find the new set of solutions. A con¬ 
tinuous set of o's results from this procedure having ’\iice" properties. 
First, the new set of solutions is well behaved in that meaningless sets 
of solutions are not encountered and convergence is quick. Second, the 
procedure gives a "visual" picture of the behavior of the solutions in 
relation to the slack variable. 

Since infinite number of modulation functions exist, some method is 
needed to select the "best" solution. The criteria used in this study 
was minimizing the sum of the squares of the magnitudes of the next set 
of current harmonics to be eliminated beyond those of present concern. 
For example, if the fifth and seventh harmonics were being eliminated 
from h( t), then the procedure is to minimize the sum of the square of 
the magnitudes of the eleventh and thirteenth harmonics. 

A flow chart illustrating optimization of the modulation function 
is given by Figure B.2. 

Figure B.3 illustrates the form of the modulation function for the 
case of elimination of fifth and seventh harmonics and for the case of 
elimination of fifth, seventh, eleventh, and thirteenth harmonics. 








L Figure B.2. Flow Chart for Optimizing 

Modulation Function 

















2.2 Modulation Function Optimization Program 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ecccecceccecceceeccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


cccccccccccccccccccccocccccceccccccccccccccccccccccccccccccccccccccccccc 
ccoceccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c n number of harmonics to be eliminated 

c n+1 harmonics are eliminated 

c r(i) - harmonics to be eliminated 

c a(i) - these are the angles from the infinite 

c search .i.e. this is the seed 

c rmax - upper limit of slack variable 

c rmin - lower limit of slack variable 

c points - number of solutions or points to be saved 

c msq - number of harmonics to be squared and sumed 

c tr(i) - harmonics to be summed and squared excluding tr(1) 

c tr(1) - any misc. harmonic that the user wishes to view 

c incre - used for increasing or decreasing the slack variable 

c incre*1 increment the slack variable from min to max 

c incresO decrement the slhck variable from max to min 

c hmag - this is the sum of sqaures 

dimension a(10),r(10),f(10),h(10),tr(10) 

common pi,conv 

call sopen(O.O) 

pi=3.141592653 

conv*180.0/pi 

read(9,*) n 

write(6, # ) ' n=' f n 

read(9,»> (r(i),i*1,5) 

write(6,*) * harmonics*',(r(i),is1,n) 

read(9,*) (a(i),1=1,5) 

wrlte(6,*) ' alphass',(a(i),is1,n) 

read(9,*) rmax 

write(6,*) ' rmax*',rmax 

read(9,*) rmin 

write(6,») * rmlns',rmln 

read(9,*) points 

wrlte(6,*) ' pointss',points 

read(9, # ) msq 

write(6,») ' msq=',msq 

read(9,*) (tr(i),i=i f 5) 

writ«(6,») ' hsquar harmonics*',(tr(i),i=2,rasq+1) 
read(9,*) incre 
write(6,») ' inerts',incre 

rlnos( rmax-rmin )/points 
maxr=int( rmax/rinc ) 
minrslnt( rmin/rinc ) 
lraaxsmaxr-minr 








do 50 i=1,5 
a(i)sa(i)/conv 
continue 

do 100 irsminr,maxr 

if( incre .eq. 1 ) r(n+1)=float(ir)»rine 
if( incre .eq. 0 ) r(n+1)sfloat( maxr+minr-ir )«rinc 
call ff( a,r,f,n+1 ) 
call square( f,fsquar,n+1 ) 
call look( a,r,f,fsquar,n+1 ) 
call shar( a,tr,h,n+1,msq+1 ) 
call shsq( h,hsquar,msq ) 
hmag=sqrt( hsquar ) 
write(1 f *) ( a(i)»conv,isl t n+1) 
write(2,1000) h(1) 
write(3,1000) h(2) 
write(4,1000) h(3) 
write(7,1000) r(n+1) 
write(8,1000) hoag 
continue 
call sclose(O.O) 
format( el0.4 ) 
format( 10(e15.7,5x) ) 
formate a=') 
format(’ harmonicas') 
stop 
end 









cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccecccccceccocccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine fa( a,da,n ) 
dimension a(10),da(10) 
common pi.conv 

do 200 i=1,n 
a(i)=a(i)-da(i) 

200 continue 
return 
end 









cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cecccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine close c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sclose(dummy) 
closed) 
close(2) 
close(3) 

closed) 
close(7) 
closed) 
close(9) 
return 
end 










cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine df c 

C C 

ccceeccccccocccceccccecccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine fdf( a,r,f,df,n ) 
dimension a(10),r(10),f(10),df(10,10) 
common pi,conv 


do 200 is1 ( n 
signs 1.0 

do 100 Jsl, n 

df(i,j)=-sign»r(i)»sin( r(i)»a(j) ) 
slgns-sign 

100 continue 

200 continue 
return 
end 
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cccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccc 
cccccecccccccccccccccccccccccccccccccccccccccccccecccccccccccccccccccccc 
c c 

c subroutine harmonics c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine shar( a,tr f h f n,m ) 
dimension a(10),h(10),tr(10) 
common pi,conv 

do 200 is1,m 
signs 1.0 
h(i)s0.0 

do 100 jsl.n 

h(i)sh(i)tJ*.0/( tr(i)*pi )*sign*cos( tr(i)»a(j) ) 
signs-sign 

100 continue 

200 continue 
return 
end 









ccccccccccccccccceccccccecceccceccccecccccccccccccccccccccccccccccoccccc 

cececccccccccccccccccccccccccccccceccecccccccccccccccccccccccccccccccccc 







ccececcccccccecccecccccccccccccccccccecccccccccccccccccccccccccccecccccc 

ceococcccceccecccccccccecceccccccccccccccceccccccccccccccccccccccccccccc 


c c 

c subroutine look c 

c c 

ecccoeccececccccccccccccccccccceccececcccccccccccccccccccecccccccccccccc 
cceeoecccccccccecccccccccccccccccccccccccccccecccccccccccccccccccccccccc 
subroutine look( a,r,f,fsquar,n ) 
dimension a(10),r(10),f(10),df(10,10),dfi(10,10) 
common pi,conv 
kount=30 
fminsO.OI 

do 2500 loopsl.kount 
call fdf( a,r,f,df,n ) 
call matv( df,n,det,dfl ) 
call mult1( f t dfi,da,n,n ) 
call fa( a,da,n ) 
call ff< a,r,f,n ) 
call square( f,fsquar,n ) 

If< loop .gt. 5 .and. fsquar .It. fmln ) go to 2900 
2500 continue 

2900 dummy=0.0 

3000 format( 10(e15.7,5x) ) 

3050 format( *0* ) 

3100 format( ' as* ) 

3150 formate ' f=' ) 

3200 formate ' fsquars* ) 

4000 formate'test',13) 

return 
end 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 


SUBROUTINE MATV 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE MATV(A,N,DET,B) 

DIMENSION INDEX(20,2),IPVOT(20),A(10,10),B(10,10),PIVOT(20) 

ONES 1.0 
ZERO=0.0 


c 

DO 10 1=1,N 

DO 5 J=1,N 
B(I,J)=A(I,J) 

5 CONTINUE 

10 CONTINUE 

DET=ONE 
DO 15 J=1,N 
15 IPVOT(J)=0 

DO 75 1=1,N 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 15 STATEMENTS C 
C FOR SEARCH FOR PIVOT ELEMENTS C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
T=ZERO 
DO 40 J=1,N 

IF(IPVOT(J)-1) 20,40,20 
20 DO 35 K=1,N 

IF(IPV0T(K)-1) 25,35,100 
25 CONTINUE 

D1=ABS(T) 

D2=ABS( B(J,K) ) 

IF(D1-D2)30,35,35 
30 IROW=J 

ICOLsK 
T=B(J,K) 

35 CONTINUE 

40 CONTINUE 
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IPVOT(ICOL)alPVOT CICOL)+1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C FOLLOWING 10 STATEMENTS PUT C 
C PIVOT ELEMENT ON DIAGONAL C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF(IROW-ICOL) 45,55,45 
45 DETx-DET 

DO 50 La1,N 
TsB(IROW,L) 

B(IROW,L)aB(ICOL,L) 

50 B(ICOL,L)sT 

55 INDEX(I,1)sIROW 

IMDEX(I,2)*IC0L 
PIVOT(I)«B(ICOL,ICOL) 

DETsDET*PIVOT(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 3 STATEMENTS TO DIVIDE PIVOT C 

C ROW BY PIVOT ELEMENT C 

C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

B(ICOL,ICOL)=ONE 
DO 60 Lxl.N 

60 B(ICOL,L)=B(ICOL,L)/PIVOT(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 7 STATEMENTS TO REDUCE C 
C NON-PIVOT ROWS C 
C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 75 LI=1,N 
IF(LI-ICOL) 65,75,65 
65 T=B(LI,ICOL) 

B(LI,ICOL)*ZERO 
DO 70 L«1,N 

70 B(LI,L)*Bai,L)-B(ICOL,L)»T 

75 CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C c 
c FOLLOWING 11 STATEMENTS TO C 
C INTERCHANGE COLUMNS C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 







$1 


DO 95 1=1,N 
LsH-I+1 

IF(INDEX(L,1)-INDEX(L,2) ) 85,95,85 
JROVfsINDEX(L, 1) 

JCOL=INDEX(L,2) 

DO 90 Ksl,N 

TsB(K.JROW) 

B(K,JROW)=B(K,JCOL) 

B(K,JCOL)=T 
CONTINUE 
CONTINUE 
RETURN 
END 









cccccccccccccccccccccccccccccccccccccccccccccccecccecccccccccccccccccccc 

cccccccccccccccccccccccccccceccccccccccccececccccceccccccccccccccccccccc 


subroutine multiply 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecccc 
ecccoccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine multi( a,b,c f iO,kO ) 
dimension a(10),b(10,10),c(10) 
do 200 isl.iO 
c(i)s0.0 
do 100 k=1,k0 

c(i)=c(i)+b(i,k) # a(k) 

100 continue 

200 continue 
return 
end 
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eceecoccccoccccccccccccccccccccccccccccccccccccccccccccccccccccccecccccc 

cccececccccoccccccccccccccccccccccccoccccccccccccccccccccccccccccccccccc 


subroutine open 


ccccccccccccccoccccoccccecccoccccccccccccccccccccccccccccccccccccccccccc 
coeececccccceeccccoocccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sopen(dummy) 
open(unit=1,file="dataO/alpha") 
rewind(1) 

open(units2,file= M dataO/fundamental") 
rewind(2) 

open(units3,files N dataO/har2 tt ) 
rewind(3) 

open(units4,file= , 'data0/har3 ,, ) 
rewind(4) 

open(units7 y flle="data0/ioghar n ) 
rewind(7) 

open(unit=8,file="data0/hmag") 
rewind(8) 

open(unit=9,file="input n ) 

rewind(9) 

return 

end 










ccoccccccccccccccccccccccccccccocccccccccccccccccccccccccccccccccccccecc 
cccccccccccccccccccccccccccccccccccccccccccccccccccceecccccecccccccccccc 
c c 

c subroutine square c 

c c 

ccccccccccccceeccccccccecccccccecccccccceccecccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccccccccccccc 
subroutine square( f,fsquar,n ) 
dimension f(10) 
fsquarsO.O 
do 200 1=1,n 

fsquar=fsquar+f(i)*f(i) 

200 continue 
return 
end 


115/116 









oecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine square c 

c c 

cccccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
occccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine square( f,fsquar,n ) 
dimension f(10) 
fsquar=0.0 
do 200 i=1,n 

fsquarsfsquar+f(i)*f(i) 

200 continue 
return 
end 
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APPENDIX C 


FOURIER SPECTRUM PROGRAM 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c frequency spectrum c 

c c 

cccccccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c tor - Input data (this is a real variable) 

c x complex variable containing the input data 

c before the fft subroutine call and 

c the complex frequency after the fft subroutine call 

c mag - magnitude of the complex frequency 

c n number of data points (must be a power of 2) 

c inv - used by fft 

c inv=0 fft 

c inv=1 inverse fft 

c maxi - used to decided highest frequency to be plotted 

c avemag - used to normalize data to 100% of the largest 

e harmonic 

dimension tor(1030) 
complex x(1030),cmplx 
real mag 
call sopen(O.O) 

pi=3.1415926535 

rms=0.0 

read(9,*) n 

write(6,*) * n=',n 

read(9,*) inv 

write(6,*) * invs'.inv 

read(9,*) maxi 

write(6,*) ' maxi=',maxi 

read(9,*) iwrite 

write(6,*) * iwrite='.iwrite 

do 100 isl.n 

if( inv .eq. 0 ) read(1,2000) tor(i) 

if( inv .eq. 0 ) x(i)=cmplx( tor(i),0.0 ) 

if( inv .eq. 1 ) read(4,2000) x(i) 

if( iwrite .eq. 1 .and. inv .eq. 0 ) write(4,») x(i) 

100 continue 

call fft( x.n.inv ) 

do 200 is1,maxi 
magscabs( x(l) ) 
phase=atan2( yy.xx ) 

avemag=amax1( cabs(x(1)),cabs(x(2)),cabs(x(3)),cabs(x(4)) ) 

rmsarms+mag ## 2 
write(3»*) mag/avemag*100.0 













cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

e subroutine close c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccceccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sclose(dummy) 
closed) 
close(2) 
close(3) 

closed) 
close(7) 
dosed) 
close(9) 
return 
end 
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'WW7TOW!trrr 


c 

o 

c 

c**The following fft subroutine is taken almost verbatim 
c**from Ahmed N., Rao K.R., Orthogonal Transforms for 
c**Digital Signal Processing, Springer-Verlag, 1975., pp79 
c 
c 
c 

subroutine fft(x,n,inv) 
c 

c this program implements the fft algorithm to compute the discrete 
c fourier coefficients of a data sequence of n points 
c 

c calling sequence from the main program: 
c call fft(x,n,inv) 
c n: number of data points 

c x: complex array containing the data sequence. In the end dft 

c coeffs. are returned in the array. Main program should 

c declare it as— complex x(1024) 

c inv: flag for inverse 

c inv=0 for forward transform 

c inv=1 for inverse transform 

c 

complex x(1030),w ( t,cmplx 
c 

c calculate the # of iterations (log. n to the base 2) 

c 

iter=0 

iremsn 

10 iremsirem/2 

if (lrem.eq.O) go to 20 
iter=iter+1 
go to 10 
20 continue 
slgns-1 

if (inv.eq.1) signal 
nxp2=n 

do 50 ltd,iter 
c 

o computation for each iteration 

c nxp:number of points in a partition 

c nxp2:nxp/2 
nxp*nxp2 
nxp2snxp/2 

wpwr*3.141592/float(nx p2) 
do 40 m*1,nxp2 









calculate the multipier 

argsfloat(m-1)»wpwr 
w=cmplx(cos(arg),sign*sin(arg)) 
do 40 mxpsnxp.n.nxp 

computation for each partition 

Jl=mxp-nxp+m 
J2=J1+nxp2 
t=x(j1)-x(J2) 
x(jl)sx(jl)+x(J2) 

40 x(J2)st»w 
50 continue 

unscramble the bit reversed dft coeffts. 

n2=n/2 

n1=n-1 

J=1 

do 65 isl.nl 
lf(i.ge.j) go to 55 
t=x(j) 
x(J)=x(i) 
x(i)=t 
55 k=n2 

60 if(k.ge.j) go to 65 
j=j-k 
k=k/2 

go to 60 
65 J=j+k 

if (lnv.eq.1) go to 75 
do 70 i»1,n 
70 x(i)sx(i)/float(n) 

75 continue 
return 
end 












cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccoccccccccc 
c c 

c subroutine open c 

c c 

ccccccccccccccocccccccccccccccccocccccccccccccccccccccccocccccccccccoccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccoccccccccc 
subroutine sopen(dummy) 
open(units1,flies"torque") 
rewlnd(l) 

open(units2,files"dataO/phase") 
rewind(2) 

open(units3,flles"dataO/Percent") 

rewind(3) 

open(units4 f files n flie1") 
rewind(4) 

open(unlts7 f fiies"data0/misc2") 
rewind(7) 

open(unlts8,files w data0/misc3") 
rewind(8) 

open(units9,files"lnput") 

rewind(9) 

return 

end 








APPENDIX D 


PERFORMANCE PROGRAM 


ccccccccccccccccccccccccccccccoccccccccccccccccceccccccocceccccccccceccc 

cccccccccccccccccccccccccccccccccccccccecccceccccccccccccccccccccccceccc 


c c 
c program to calculate developed power- c 
c c 
c speed points for brushless dc motor c 
c c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccc 

c 

c 

c rpm - speed of dc motor 

c wm - angular speed of dc motor 

c p number of poles 

c t time 

c kO number of intervals between print outs 

c h interval of integration 

c w electrical angular frequency of pmm 

c ws electrical angular frequency of pmg 

c alpha - commutation advance angle (degrees) 

c ra stator resistance 

c xla - stator leakage Inductance 

c rO - choke coil resistance 

c xlo - choke coil inductance 



xk - motor constant 

xkg - generator constant 

vd - inverter input voltage 

x - array of stator phase currents 

tp - period 

ang - electrical angular displacement 

xlamda - electrical angular displacement (pmg) 

m - number of phases 

pd - total developed power 

td - total developed torque 

pdl - phase 1 developed power 

pd2 - phase 2 developed power 

pd3 - phase 3 developed power 

e(1) - excitation voltage of phase 1 

e(2) - excitation voltage of phase 2 

e(3) - excitation voltage of phase 3 


dimension f(6),v(6),x(6),xx(6),a(6),r(6),e(6),alphs(2) 

common rr,rd,r0,ra,r(6),pi,xla,xl0,e(6),vm,vdc,ivdon.imodon 

real lav.irms 

dummy=0.0 

call popen(dummy) 

pi*3.1415926536 

conv=pi/l80.0 

read(9,*) lstop 
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write(6,*) * istops•,istop 

read(9,*) alpha 

write(6,*) * alpha**,alpha 

read(9,») ( a(i),isl,5 ) 

write(6,*) • as*,( a(i),isi,5 ) 

read(9»*) h 

writer, 11 ) ' hs'.h 

read(9,*) rO 

write(6,») * rOs',rO 

read(9,») xlO 

wrlte(6,*) • xlOs*,xlO 

read(9,*) ra 

write(6,») * ras’.ra 

read(9,*) xla 

write<6,») » xlas'.xla 

read(9,») rd 

write(6,*) » rdr'.rd 

read(9,*) rpmO 

write(6, # ) * rpmO=* t rpmO 

read(9,*) xk 

write(6,») • xk=’,xk 

read(9,*) xkg 

write(6,*) ' xkgs'.xkg 

read(9»*> pmax 

write(6,*) ' pmaxs*,pmax 

read(9 t *) pmin 

write(6,*) * pmins'.pnin 

read(9,*) vdc 

write(6,*) ’ vdc=',vdc 

read(9,*) points 

write(6,*) ' points**.points 

read(9,*) ivdon 

write(6,•) » ivdons*,ivdon 

read(9 # *) imodon 

write(6,*) » imodons*.imodon 

read(9,*) ichon 

write(6,») ' ichons*.ichon 

read(9,*) iwgon 

write(6,*) ’ iwgonsiwgon 

read(9,») ( x(i),i*1,3 ) 

write(6,») » xs»,( x(i),isi,3 ) 

read(9,*) iarea 

write(6,*) * iareasiarea 

read(9,*) xwg 

write(6,*) ' xwgs'.xwg 

read(9,*) ifull 

write(6,») • ifulls*,ifull 

do 4 J*1,5 

if( imodon .eq. 1 ) a(J)s a(J)»pi/180.0 
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1 


0 ) a(j)spi/2.0 

0 ) a(1)*30.00*pi/180.0 


) wg=7950*pl 
) wg=xwg*wm 


if( imodon .eq. 
lf( imodon .eq. 
continue 
m=3 

pfw=CrpmO/45000.0)«»2«330.0 
wm=rpm0*pi/30.0 
if( lwgon .eq. 1 
if( iwgon .eq. 0 
umsabs(wm) 
p=4.0 
t=0.0 

w=p/2.0*wm 

ws=2.0*wg 

vm=sqrt(3.0)*xkg*wg 
write(6,*) * vms'.vm 
alpha=alpha*conv 

area=0.5+0.5*( pi/3.0-a(1)+a(2)-a(3)+a(4)-a(5) )/( pi/3.0 ) 


0 .and. iarea .eq. 1 ) vdcsvdc/area 
eq. 1 ) alphs(1)=acos( (vdc*aqrt(2.0))/(1.0*vm*1.35) ) 
eq. 1 ) alpha(2)=acoa( (vdc*aqrt(2.0))/(2.0*vm*1.35) ) 
eq. 0 ) alpha(1)*acoa( (vdc*aqrt(2.0))/(0.5*vm*1.35) ) 
eq. 0 ) alpha(2)sacos( (vde*aqrt(2.0))/(1.0*vm*1.35) ) 

* alpha(1)=' ,alphad)/conv 

* alphs(2)=‘,alpha(2)/conv 
' area=',area 

* f*',w/(2.0*pi) 

' fa=' ,wa/(2.0*pi) 

' mag of els',xk*w/2.0 


ifC imodon 
if( ifull 
if< ifull 
if( ifull 
if( ifull 
write(6,*) 
write(6,*) 
write(6,*) 
write(6,•) 
write(6,*) 
write(6,*) 
xxd )=0.0 
xx(2)=0.0 
xx(3)=0.0 
pdsO.O 
pinsO.O 
iavsO.O 
irmasO.O 
e1ave=0.0 
vdavesO.O 
eIrmasO.O 
sumsO.O 
tp*2.0*pi/w 
tminspain*tp 
tmaxspmax*tp 
tmxmnstmax-tmin 
soale*30.0/pointa 

kO»int( tmxmn/tp*acale/aba(rpmO)/h ) 

kountsO 

index*1 

rewind(9) 

read(9, # ) latop 

lf( latop .eq. 1 ) go to 305 







kountskount+1 

ang=w*t 

iprinl=int(ang/2.O/pi) 

angsang-floatUprinl )*2.0»pi 

e(1)=xk*w § 2.0/p*sin(ang-alpha) 

e(2)sxk»w*2.0/p*sin(ang-alpha-2.0/3.0«pi) 

e(3)=xk # w*2.0/p»sin(ang-alpha+2.0/3.0»pi) 

wstsws*t0 

if( wst .ge. 2.0*pi ) tOsO.O 
WStswa*tO 

if( ifull .eq. 0 ) call half( wst.xlaroda ) 
if( ifull .eq. 1 ) call full( wst.xlamda ) 
call fv( a,ang,v,vd,alphs,xlamda ) 
call ff( index,t,f,x,xi,ang,del,v ) 
if( ichon .eq. 1 ) call check( m,x,f ) 
call derk( ra,x,f,t,t0,h,index ) 
if( index .eq. 2 ) go to 200 
if( t .It. tmin ) go to 100 
f(1)s( x(1)-xx(1) )/h 
fC2)=( x(2)-xx(2) )/h 
f(3)=C x(3)-xx(3) )/h 
do 86 kk=1,3 
xx(kk)=x(kk) 
continue 

call ffi( fi,f,ang ) 

call fvap( wst,vap,alphas ) 

v1p=vd-fi*xl0/2.0-xi*r0-x(1)»(r(1>-ra) 

vscr=vap-v1p 

if( ang .It. 3.0/3.0»pi .and. wst .It. 2.0/3.0»pi ) 
vcsr=x(1)»(r(1)-ra) 

if( ang .ge. 3.0/3.0»pi .and. wst .ge. 3.0/3.0»pi .and. 

wst .It. 4.0/3.0»pi > vscr=x(1)*(r(1)-ra) 
ifC t .ge.tmax ) go to 300 
if( kount .It. kO ) go to 100 
suassum+l.0 
pd1sx(1)*e(1) 
pd2=x(2)*e(2) 

pd3*x(3)*e(3) 

td=( pd1+pd2+pd3 )/wm 

pd=pd+( pd1+pd2+pd3 ) 

pin*pin+( v(1)»x(1)+v(2)«x(2)+v(3) § x(3) ) 

irms=irms+x(1) # *2 

iav«iav+abs( x(1) ) 

e 1 ms*e 1 rms+e (1) ••2 

e1ave*e1ave+abs( e(1) ) 

vdave*vdave*vd 

hsin6*hsin6+td*sin(6.0*ang) 

hcos6*hcos6+td*cos(6.0*ang) 

hsin12*hsin12+td # sin(12.0»ang) 

hcoa12»hcos12+td*cos(12.0*ang) 







hain18=h3in18+td*ain(18.0*ang) 

hcoa18=hcos18+td*cos(I8.0*ang) 

hsin24shain24+td*sin(24.0»ang) 

hcoa24=hcos24+td*coa(24.0*ang) 

txst-tmxmn 

txxst*w»180.0/pi 

write(1,2000) x(1) 

wrlte(2,3000) txx 

write(3,2000) td 

write(4,2000) v(1) 

write(7,2000) vd 

write(8,2000) x(1)+x(2)+x(3) 

kountsO 

lf( t .It. tmax ) go to 100 
305 pdspd/aum 
pinspin/aun 
Irmasirma/sum 
irma=sqrt( Irma ) 
iavslav/sun 
e1rma=e1roa/aun 
elrmasaqrtC elrma ) 
elaveselave/aum 
vdave=vdave/aum 
haln6=hsln6/sum 
hcos6=hcos6/sum 
h6=hsin6**2+hcos6**2 
h6=aqrt( h6 )/6.0 
hain12shain12/sum 
hcoa12=hcos12/sum 
h12shaln12**2+hcoa12**2 
h12=aqrt( h12 )/12.0 
hslnl8shsln18/sum 
hcoa18=hcos18/sum 
hl8=hslnl8**2+heoa18**2 
h18=aqrt( h18 )/18.0 
hsln24shaln24/sum 
hcos24shcoa24/sum 
h24ahsin24*«2+hcos24**2 
h24ssqrt( h24 )/24.0 
600 taavs(pd-pfw)/wm 

wrlte(6,*) ' taavs',tsav 
effs(pd-pfw)/pin*100.0 
lf< eff .It. 100.0 ) go to 710 
effspin/(pd+pfw)*100.0 
710 wrlte(6,*) ' pins',pin 

wrlte(6,*) ' effs',eff 
write(6,*) ' lavs',lav 
write(6,*) ' Irmas',irms 
wrlte(6,*) ' elave*',e1ave 
write(6,*) ' vdaves',vdave 
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write(6,») ' h6=',h6 
write(6,») • h12s',h12 
write(6,») • h18s’,h18 
write(6,») • h24s*,h24 
1000 format( 5(e12.4,5x) ) 

2000 formate e10.4 ) 

3000 formate f8.3 ) 

call pclose(dummy) 

stop 

end 
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cccccecocccocecceccccecccccccccccccecccccccccccccccccccccccccccccccccccc 
ccccccccccceccccccccccccccccccccccccccccccccccccccccccccccccccccccceccec 
e c 

c subroutine check c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine check( m,x,f ) 
dimension x(6),f(6) 

common rr,rd,rO f ra,r(6),pi,xla,xlO,e(6),vm,vdc,lvdon t imodon 
do 10 i=1,m-1 

slgnf=abs( f(l) )/f(l) 
signx=abs( x(l) )/x(i) 

if( abs( f(i) ) .gt. 1.0e08 ) f(i)=signf»1.0e08 
ifC abs( x(i) ) .gt. 1.0e06 ) x(l)=signx # 1.0e06 
10 continue 

f(3)=-( f(1)+f(2) ) 

return 

end 






cceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CCCCCCCCCCCCCOCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


i 


c subroutine pclose c 

c c 

cccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine pclose(dummy) 
close(1) 
close(2) 
close(3) 
close(4) 
close(7) 
dose(8) 
close(9) 
return 
end 









cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

c subroutine derk c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c derk is a fourth-order, fixed increment runge-kutta 

c integration routine 

c 

c m 

C X 

c f 

c h 

c t 

c index 

c 1. if 

c 2. if 

c 

subroutine derk( m,x,f,t,tO,h,index ) 
dimension x(6),f(6),q(400) 
if( index .eq. 2 .) go to 19 

18 kxxsO 
index=2 
do 35 i=1,m 
j=i+300 

35 q(j)=x(i) 

19 kxxakxx+1 
go to (1,2,3*4),kxx 

1 do 5 i*1,m 
q(i)=f(i)*h 

5 x(i)=x(i)+q(i)/2.00 
tat+h/2.00 
tOstO+h/2.00 
return 

2 do 6 isl.m 
Jai+100 
kai+300 
q(J)af(i)»h 

6 x(i)aq(k)+q(J)/2.00 
return 

3 do 7 ia1,m 
Jai+200 
kai+300 
q(J)af(i)»h 

7 x(i)aq(k)+q(J) 
tat+h/2.00 
t0at0+h/2.00 
return 


- number of simultaneous differential equations 

- array of dependent variables 

- array of derivatives of independent varibles 

- increment 

- time 

- indicator 

exit routine with indexal, solution found at=t+h 
exit routine with indexa2, go back to re-evaulate 





do 8 i=1,m 
jsi+100 
kai+200 
1 * 1+300 

x(i)=q(l)+( ( q(i)+2.00*q(k)+2.00*q(k)+f(i)*h )/6.00 ) 

indexsl 

return 

end 
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cccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


cccccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine ff( lndex,t,f f x,xl,ang,del,v ) 
dimension f(6),v(6),x(6) 

common rr,rd,rO f ra,r(6),pi,xla,xl0,e(6),vm,vdc.ivdon,imodon 
rrsrd+ra 

if( index .eq. 1 ) r(1)=rr 
lf( Index .eq. 1 ) r(2)=rr 
lf( index .eq. 1 ) r(3)=rr 

if( ang .ge. 0.0 .and. ang .It. 1.0/3.0*pi ) goto 10 

if( ang .ge. 1.0/3.0»pi .and. ang .It. 2.0/3.0»pi ) goto 20 

if( ang .ge. 2.0/3.0«pi .and. ang .It. 3.0/3.0»pi ) goto 30 

if( ang .ge. 3.0/3.0»pi .and. ang .It. 4.0/3.0»pi ) goto 40 

if( ang .ge. 4.0/3.0«pi .and. ang .It. 5.0/3.0»pi ) goto 50 

if( ang .ge. 5.0/3.0«pi .and. ang .It. 6.0/3.0»pi ) goto 60 

10 n1=1 

n2=3 
n3=2 

if( icycle .eq. 1 ) iflag=0 
itests10 
icyclesO 
go to 100 
20 nisi 

n2=2 

n3»3 

iteat=20 
go to 200 
30 n1s2 

n2si 

n3*3 
itest*30 
go to 100 
10 n1s2 

n2*3 

n3*1 

itests40 
go to 200 
>0 n1s3 

n2*2 
n3*1 

iteat«50 
go to 100 









n2=1 
n3=2 
itests60 
lcyclesl 
go to 200 

100 lf( x(n2) .It. -0.001 .and. index .eq. 1 ) iflag=itest 

if( x(nl) .It. -0.001 .and. index .eq. 1 ) r(n1)srr+1000.0*abs( x(nl) ) 

if( x(n3) .gt. 0.001 .and. index .eq. 1 ) r(n3)*rr+1000.0*abs( x(n3) ) 

if( iflag .eq. itest .and. index .eq. 1 ) r(n2)=rr+1000.0*abs( x(n2) ) 

do 150 is1,3 

if( r(i) .gt. 1000.0 ) r(i)=1000.0 
150 continue 

xlsxla 

f(n1)s( -x(n1)*r(n1)+v(n1)-e(n1) )/xl 
f(n2)=( -x(n2)*r(n2)+v(n2)-e(n2) )/xl 
f(n3)=( -X(n3)*r(n3)+v(n3)-e(n3) )/xl 
return 

200 if( x(n2) .gt. 0.001 .and. index .eq. 1 ) iflag=itest 

lf( x(n3) .gt. 0.001 .and. index .eq. 1 ) r(n3)=rr+1000.0»abs( x(n3) ) 

if( x(nl) .It. -0.001 .and. index .eq. 1 ) r(n1)xrr+1000.0«abs( xCnl) ) 

if( iflag .eq. itest .and. index .eq. 1 ) r(n2)=rr+1000.0*abs( x(n2) ) 

do 250 i=1,3 

if( r(i) .gt. 1000.0 ) r(i)=1000.0 
250 continue 

xl=xla 

f(n1)=( -X<n1)*r(n1)+v(n1)-e(n1) )/xl 
f(n2)s( -x(n2)*r(n2)+v(n2)-e(n2) )/xl 
f(n3)*( -x(n3)*r(n3)+v(n3)-e(n3) )/xl 

return 

end 







cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccceecceccccccecccccccccccccccccccccccccccccccccccccccccccccccecccccccc 
c c 

o subroutine fi c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccoeececcccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine ffi( fi,f,ang ) 
dimension f(6) 

common rr,rd,r0,ra # r(6),pi t xla,xlO,e(6),vm,vdc,ivdon,imodon 
if( ang .ge. 0.0/3.0«pi .and. ang .It. 1.0/3.0*pi ) fi=-f(2) 

if( ang .ge. 1.0/3.0»pi .and. ang .It. 2.0/3.0*pi ) fi* f(1) 

if( ang .ge. 2.0/3.0»pi .and. ang .It. 3.0/3.0«pi ) fi=-f(3) 

if( ang .ge. 3.0/3.0»pi .and. ang .It. 4.0/3.0*pi ) fi* f(2) 

if( ang .ge. 4.0/3.0«pi .and. ang .It. 5.0/3.0*pi ) fi=-f(1) 

if< ang .ge. 5.0/3.0»pi .and. ang .It. 6.0/3.0*pi ) fis f(3) 

return 
end 








imu. i i mm u i wm p. 


cccccccccceccccccocccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccc 
c c 

c subroutine full-bridge c 

c c 

cccccccceccccccccccceecccccccoccccecoccccccccccccccocccccccccccccccccccc 
cecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine full( wst,xlamda ) 

common rr,rd,r0 t ra,r(6),pi,xla,xl0,e(6),vm,vdc,ivdon f imodon 


if( wst .ge. 0.0 
wst-0.0 

ifC wst .ge. 1.0»pi/3.0 
wst-1.0*pi/3.0 
if( wst .ge. 2.0»pi/3.0 
wst-2.0*pi/3.0 
ifC wst .ge. 3.0»pi/3.0 
wst-3.0*pi/3.0 
if( wst .ge. 4.0*pi/3.0 
wst-4.0»pi/3.0 
if( wst .ge. 5.0»pi/3.0 
wst-5.0*pi/3.0 
return 
end 


.and. wst .It. 1.0*pi/3.0 ) xlamda= 
and. wst .It. 2.0*pi/3.0 ) xlamda= 
and. wst .It. 3.0*pi/3.0 ) xlamdas 
and. wst .It. 4.0*pi/3.0 ) xlamda= 
and. wst .It. 5.0*pi/3.0 ) xlamda= 
and. wst .It. 6.0*pi/3.0 ) xlamdas 







cccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccc 
cceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 
c subroutine half-bridge c 
c c 
cccccccccccccceececccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccecccccccccccccccccccecceccccccccccccccccccccccccccccccccccccccc 


subroutine half( wst.xlamda ) 

common rr,rd,rO,ra,r(6),pi,xla,xlO,e(6),vm,vdc,ivdon,imodon 
if( wst .ge. 0.0 .and. wst .It. 2.0*pi/3.0 ) xlamda= 

c wst-0.0 

if( wst .ge. 2.0*pi/3.0 .and. wst .It. 4.0*pi/3.0 ) xlamda= 
c wst-2.0»pi/3.0 

if( wst .ge. 4.0*pi/3.0 .and. wst .It. 6.0»pi/3.0 ) xlamda= 
c wst-4.0*pi/3.0 
return 
end 
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open(unit=3,file= B dataO/Torque") 
rewind(3) 

open(unit=4,file="dataO/Van") 
rewind(4) 

open(unit=7,file="dataO/Vd") 
rewind(7) 

open(unit=8,file="dataO/In"> 
rewind(8) 

open(unit=9,file*"input") 

rewind(9) 

return 

end 
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cceeccccccccceecceeccoccccccccecccccccceeccccccccccccccccccccccccccccccc 
ccccoccoceeccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine vd c 

c c 

cceeccceocceoccccecececceccccccoecceeccccccccccccccccccccccecccccccccccc 
ccccccccccececcccccccccccccocccccccccccccccccccccccccccccccccccccccccccc 
subroutine fv( a,ang t v,vd,alphs,xlamda ) 
dimension a(6),v(6),alphs(2) 

common rr,rd,r0,ra,r(6),pi,xla,xl0,e(6),vm,vdc,ivdon,imodon 
beta=xlamda+alphs(1) 

if( ivdon .eq. 1 .and. ifull ,eq. 1 ) vd=vm*sin( beta+pi/3.0 > 

if( ivdon .eq. 1 .and. ifull .eq. 0 ) vd=vm*sin( beta+pi/6.0 ) 

if( ivdon .eq. 0 ) vd=vdc 

shift=30.0/180.0»pi 

theta=shift+ang 

call svxn( vx,theta,a ) 

van=vx 

theta=shift+ang-2.0/3.0 # pi 
call svxn( vx,theta,a ) 
vbnsvx 

theta*shift+ang+2.0/3.0«pi 
oall svxn( vx,theta,a ) 
vcnavx 

v(1)*vd*van 

v(2)avd»vbn 

v(3)*vd»vcn 

return 

end 







CJVkUUAH WP. f JUWW JM JLP.Y A 1 A AA \ 


w* ..' ‘ K - i. m K m V 



subroutine fvap( wst,vap,slps ) 

common rr,rd,rO,ra,r(6),pl,xla,xlO,e(6),vm,vdc,ivdon,imodon 
vaasO.O 


vab* vm*sln( wst+1.0/3.0»pi+alpa ) 
vac=-vm*sln( wst+5.0/3.0*pl+alps ) 


if( 

wst 

.ge. 

0.0/3.0»pi 

.and. 

wst 

.It. 

1.0/3.0»pi 

) vap* 

vab 

if( 

wst 

.ge. 

1.0/3.0«pi 

.and. 

wst 

.It. 

2.0/3.0*pl 

) vap* 

vac 

if( 

wst 

.ge. 

2.0/3.0»pi 

.and. 

wst 

.It. 

3.0/3.0*pi 

) vap* 

vac 

if( 

wst 

.ge. 

3.0/3.0»pi 

.and. 

wst 

.It. 

4.0/3.0»pi 

) vap* 

vaa 

if( 

wst 

.ge. 

4.0/3.0»pi 

.and. 

wst 

.It. 

5.0/3.0«pl 

) vap* 

vaa 

if( 

wst 

.ge. 

5.0/3.0*pi 

.and. 

wst 

.It. 

6.0/3.0»pl 

) vap* 

vab 


return 

end 
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ccccccccccccccccececcccecccccceecccccccccccccccccccccccccccccccccccccccc 

ccccccecccccccccccccccccccccccccccccccccccccccccccccccccceccccccccccccec 


c subroutine vxn 

c 


c 

c 

c 


ccccccccoccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine svxn( vx,th,a ) 
dimension a(6) 

common rr,rd,r0 # ra,r(6),pi,xla,xl0,e(6),vm,vdc,ivdon,imodon 

vxsO.O 

xa=0.0 

xbspi 

xcs2.0*pi 

a(6)=pi/2.0 

th=th+2.0*pi 

th=th-float( intCth/xc) )»xc 

if( th .ge. xa+a(1) .and. th .It. xa+a(2) ) vx=1.0 

if( th .ge. xa+a(3) .and. th .It. xa+a(4) ) vx=1.0 

ifC th .ge. xa+a(5) .and. th .It. xa+a(6) ) vx=1.0 

if( th .ge. xb-a(6) .and. th .It. xb-a(5) ) vx=1.0 

if( th .ge. xb-a(4) .and. th .It. xb-a(3) ) vx=1.0 

if( th .ge. xb-a(2) .and. th .It. xb-a(1) ) vx=1.0 

if( th .ge. xb+a(1) .and. th .It. xb-<-a(2) ) vx=-1.0 

if( th .ge. xb+a(3) .and. th .It. xb+a(4) ) vx=-1.0 

if( th .ge. xb+a(5) .and. th .It. xb+a(6) ) vx=-1.0 

if( th .ge. xc-a(6) .and. th .It. xc-a(5) ) vx=-1.0 

ifC th .ge. xc-a(4) .and. th .It. xc-a(3) ) vx=-1.0 

if( th .ge. xc-a(2) .and. th .It. xc-a(1) ) vx=-1.0 


return 

end 







